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The existence of more than one steady-state in a many-body quantum system driven out-of- 
equilibrium has been a matter of debate, both in the context of simple impurity models and in the 
case of inelastic tunneling channels. In this Letter, we combine a reduced density matrix formalism 
with the multilayer multiconfiguration time-dependent Hartree method to address this problem. 
This allows to obtain a converged numerical solution of the nonequilibrium dynamics. Considering 
a generic model for quantum transport through a quantum dot with electron-phonon interaction, 
we prove that a unique steady-state exists regardless of the initial electronic preparation of the 
quantum dot consistent with the converged numerical results. However, a bistability can be observed 
for different initial phonon preparations. The effects of the phonon frequency and strength of the 
electron-phonon couplings on the relaxation to steady-state and on the emergence of bistability is 
discussed. 



The existence of a unique steady-state in strongly cor- 
related quantum systems out-of-cquilibrium has been a 
subject of great interest and controversy. For the case 
of the Anderson impurity model, it has been argued us- 
ing the Bethe ansatz that a single steady-state solution 
exists [lj. However, recent calculations of the nonequilib- 
rium current based on time-dependent density functional 
theory seem to indicate that at long times the system 
reaches a dynamical state characterized by correlation- 
induced current oscillations [2 . Similarly, questions re- 
garding hysteresis, bistability and the dependence of 
the steady-state current on the initial occupation have 
been raised in the context of inelastic transport through 
nanoscale quantum dots jSHO]. 

Addressing the issue of a unique steady-state is a chal- 
lenging task for theory, as systems exhibiting bistabil- 
ity involve strong electron-electron or electron-phonon 
correlations. Under these conditions, an exact solu- 
tion is unavailable, and one has to resort to approx- 
imate methods or to numerical techniques. The for- 
mer are based on either a mean-field approximation 
or on a perturbative scheme, where the inclusion of 
higher order corrections is not always clear or system- 
atic, and thus may lead to questionable results. Numeri- 
cally brute-force approaches, such as time-dependent nu- 
merical renormalization-group techniques |10H12| . itera- 
tive |131fl"5] or stochastic |16ffl9] diagrammatic methods, 
and wave function based approaches [20J , have been very 
fruitful, but are limited in the range of parameters and 
timescales that can be studied. 

In this Letter, we address the problem of a unique 
steady-state in a system with electron-phonon couplings. 
To study the nonequilibrium dynamics, we develop an 
approach based on a reduced density matrix (RDM) for- 



malism, which requires as input a short-lived memory 
kernel |21) . We show that if a steady-state exists then 
it must be unique, regardless of the initial electronic 
preparation. However, a bistability can develop for dif- 
ferent initial phonon states. The relaxation to steady- 
state and the appearance of the bistability depends on 
the phonon frequency and the strength of the electron- 
phonon couplings. To illustrate this, we combine the 
approach with the multilayer multiconfiguration time- 
dependent Hartree method (ML-MCTDH) [H] to numer- 
ically converge the memory kernel at short times until it 
decays, and infer from it the dynamics of the system at 
all times and the approach to steady-state. 

We consider a generic model for charge transport 
through a quantum dot with electron-phonon interaction. 
The model is described by the Hamiltonian H = Hs + 
Hb + Vsb, where Hs = £d<$d is the system (quantum 
dot) Hamiltonian with creation/annihilation fermionic 
operators d^/d and energy Ed, Hb = H + H „jj where 

Hi = X^fceL n £ k a \ a k represents the noninteracting leads 
Hamiltonian with fermionic creation/annihilation opera- 
tors a\/a k , and Hp h = (b\jb a + |) represents the 

a 

phonon bath with creation/annihilation bosonic opera- 
tors b' a /b a for phonon mode a with energy u> a . The 
coupling between the system and the baths is given by 
V S b = V + V ph where Vi = £ (t k da\ + t%a k dA is the 

k£L,R V ' 

coupling between the system and the leads with couplings 
strength tf., and Vgh = d^d^2,M a (\y a + 6 a ) is the cou- 

a 

plings between the system and the phonon bath, where 
M a is the electron-phonon couplings to mode a. 

The coupling strengths were determined by various 



2 



spectral functions. The dot-leads coupling terms were de 
tcrmined from a spectral function of the form Tl r(s) = 
^E k eL, R \tk\ 2 5(e - s k 

was employed: Tl ji (e) 
a = 0.2eV and 6 =' leV. 



where a tight-binding form 



(e - Ml.-r) 2 witn 
fJ-L.R is the chemical potential 
of the left (L) or right (R) lead, respectively. Similarly, 
the electron-phonon couplings were determined from a 
phonon spectral function J(uj) — § X) Q TT'^i'^ _ 
where J (w) = |f/ue~^ is taken to be of Ohmic form. 
The dimensionless Kondo parameter, ij — deter- 
mines the overall strength of the electron-phonon cou- 
plings where A = Yl ^f- = - J — J{u>) is the reorganiza- 
tion energy (or polaron shift) and uj c is the characteristic 
phonon bath frequency. 

Following the derivation outlined in Ref. |2D for the 
Anderson impurity model, the equation of motion for the 
RDM, a(t) = Tr B {p(t)}, is given by 

ih~a (t) = C s a (t) +#(t)-~J dm (r) a (t - r) (1) 

where £$ = [Hs,---] is the system's Liouvillian, 
Tr B {- • • } is a trace over the baths degrees of freedom 
(leads and phonon baths) and p(t) is the full density ma- 
trix. In the above, 



0(t) = Tr B {iCve-* QCt Qp(0)} 



(2) 



depends on the choice of initial conditions and vanishes 
for an uncorrelated initial state (which is the case dis- 
cussed below), i.e. when p(0) = er(0) <g> p B {0), where 
<t(0) and /9b (0) are the system and baths initial density 
matrices, respectively, and C v = [Vsb, • • • ]■ The memory 
kernel, which describes the non-Markovian dependency 
of the time propagation of the system, is given by 



n(t) = Tr B {c v e-* QCt Q£p B } 



(3) 



where Q = 1 - P, P = p B (0)Tr B {- ■ ■ } and £ = [H, ■ ■ ■ ]. 

To obtain a(t), one requires as input the super-matrix 
of the memory kernel, which can be expressed in terms 
of a Volterra equation of the second type, removing the 
complexity of the projected dynamics of Eq. |3| : 

i r* 

K (t) = ih® (t) -<5>{t)£ s + - / dr$ (t - r) k (t) (4) 



with 



$(t) =Tr B {£ye-s £t ps} 



(5) 



Evaluating the super-matrix $(i) requires a tedious cal- 
culation similar to that carried out in Ref. [2D The diago- 
nal matrix elements $>u, m m(t) are given by ^3 {ip mm (t)} , 
where 



Pb (n\ ^2 

k 



t k d(t)a{(t)\m)\. (6) 



The indices i, j, m, and n can take the values 1 or 0, cor- 
responding to an occupied or an unoccupied dot, respec- 
tively. The off-diagonal elements of $ij ynm (t) are given 
by 6 jl 5 i0 ilJ mn (t) +<5,- <Wnm(*)> where 



ipmn{t) = Tr B < p B (n\ y^t k al(t) 

{ k 

- dHt)^M a (bl(t) + b a (t)) 



(7) 



Since both the operator d(t)a k (t) appearing in the 
equation for tp m n{t) anc > the full Hamiltonian conserve 
the total particle number, y> mn (i) is nonzero only for 
m = n, in which case it has a simple physical interpreta- 
tion as the time derivative of the dot population. ip mm {t) 
can be expressed in terms of the sum of the left and right 
currents eip mm (t) = I^(t) + I^,(t), where: 



rL,R 



(t) 



2e ( 
h 



keL,R 



t k {m\d(t)a\{t)\m), (8) 



and e is the electron charge. In addition, one can show 
that ipmn(t) is nonzero only for m ^ n. As a result, the 
populations and coherences of a(t) are decoupled |21| . 
and if one is interested in the behavior of the populations 
alone, only tp mm (t) is required. The resulting equations 
of motion for the diagonal elements of a(t) are given by: 



d , , 



drKi 



(t) cr mm (t-r). (9) 



If a(t) has a steady-state solution as t — > oo, then 

§i a H (*) = °- and S /o°° dTK U,mm (t) cr mm (t-r) = 

m 

can be replaced by a linear set of equations given by: 



o. 



(10) 



J °° dTK liymm (r) and cr. 



where /C,„ 

oo). The matrix JC in Eq. (10) has two eigenvalues. The 
first is Ao = corresponding to ooo = K^+kn anc ^ 
1 — cm = v — ■ For a physical steady-state 



'Coo+'Ci 



solution, both /Coo an d must share the same sign. 
Otherwise, the diagonal elements of a cannot both be 
positive. The other eigenvalue is Ai = /Coo + /Cii corre- 
sponding to coo — — en — 1- which can never represent 
a physical solution. Furthermore, since the steady-state 
depends only on the value of /Coo an d ICu a nd since both 
are independent of the initial dot population, the steady- 
state is independent of the initial preparation of the dot 
occupation and therefore, is unique. 

The steady-state solution is unique with respect to the 
electronic initial preparation. This is a strong statement 
by itself, but it does not rule out bistability for different 
initial phonon preparations. To address this question, 
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we combined the formalism described above for the RDM 
with the ML-MCTDH approach in second quantized rep- 
resentation (SQR) p0lB5]. The ML-MCTDH-SQR pro- 
vides a tool to compute the currents in Eq. ([8| numerically 
exactly. The kernel n(t) is then obtained by numerically 
solving Eq. Q. In comparison, for most model param- 
eters studied in this work, it is practically impossible to 
obtain converged values for the RDM directly from the 
ML-MCTDH-SQR, since the time to reach a steady-state 
solution is significantly longer than the maximum simu- 
lation time reachable by the ML-MCTDH-SQR. How- 
ever, since the memory kernel decays on much shorter 
timescales compared to the RDM itself [21], it is rather 
straightforward to calculate it using the ML-MCTDH- 
SQR and then solve Eq. @ for the RDM. 

To characterize the population dynamics, 
we start with a factorized initial condition of 
the form p (0) = a (0) ® p ph (0) <g) Pi ea d s (0), 
where a (0) determines whether the electronic 
level is initially occupied/unoccupied, p^ (0) = 



-p f 5> Q (bib a + 1) + J2$ a Q& + b a ) 



rep- 



initial density matrix of the phonon 



exp 



-p 



exp 

resents the 

bath. Hereby two values of the parameters 8 a are 
considered: S a — (corresponding to a phonon 
initial state equilibrated with an unoccupied dot) 
and S a — \/2uj a \ (corresponding to phonons 
equilibrated to an occupied dot). Pleads (^) = 

(sk - a\a k + ( £ k - Mi?) a \ a k 

\keL keR 

determines the initial density matrix for the leads. In the 
above P = is the inverse temperature. In all results 
shown below we take T = and \xl = —^r = 0.05eV. 

In Fig. [I] we plot the time evolution of <7n(i) (lower 
panels) and the corresponding nonzero elements of the 
memory kernel (upper panels), for two different initial 
vibrational preparations. We show the time evolution of 
<7n(f) for different values of the cutoff time t c at which 
we assume that the memory kernel has essentially de- 
cayed to zero, such that it can be safely truncated. For 
5 a = 0, it is safe to truncate the memory kernel at 
t c > 30fs while S a = \/2u a \ requires a larger cutoff time 
of t c > 80fs. Comparing the time it takes for the mem- 
ory kernel to decay (upper panels of Fig. [T]) with the time 
taken by the RDM to reach steady-state (corresponding 
lower panels of Fig. [IJ, it is clear that the latter is slower 
by nearly an order of magnitude and in some cases even 
more. Since the calculation of the memory kernel using 
the ML-MCTDH-SQR is by far the most time consum- 
ing portion of the calculation, the combination with the 
RDM formalism provides a significant saving, and more 
importantly extends the ML-MCTDH-SQR approach to 
regimes inaccessible by direct application. 

The inset of Fig. [TJshows the steady-state value of <j\\ 
as a function of 1 jt c for S a = \ / 2uj a X. For large values of 
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Figure 1: The occupation of the quantum dot <7n(i) (lower 
panels) for different cutoff times and the nonzero elements of 
the memory kernel (upper panels) for 8 a = (left panels) and 
S a — \/2oj a A (right panels). The inset shows the steady-state 
value of an for the case of S a — ^J2uj a A as a function of l/t c . 
The model parameters used are: Ed = 0.5eV, lu c — 500cm -1 , 
and A = 3000cm -1 . 



l/t c (short cutoff times) we find that the formalism may 
lead to unphysical situations in which an becomes neg- 
ative. Of course, this is expected, since only when the 
memory kernel has decayed to zero does the cutoff ap- 
proximation provide a physical meaningful solution. As 
l/t c decreases an converges and approaches a plateau 
value. In the present case of parameters, the steady- 
state value of (Tn computed for the two initial vibra- 
tional states roughly coincides. However, the dynamics 
and timescales to relax to the steady-state are clearly 
sensitive to the initial vibrational preparation. 

In Fig. [2] we plot a\\{t) for four different values of ui c 
and compare the time dependence for four different initial 
conditions, corresponding to an initial empty (n^ = 0) or 
occupied (n<j = 1) dot and to S a — or 5 a = \/2u) a A. For 
a given S a , we find that a nit) has a unique steady-state 
solution regardless of the initial dot occupation. This 
numerically converged result is consistent with the ana- 
lytical proof given above. In contrast, for two different 
initial vibrational states, a clear bistability is observed 
and the RDM decays to a different steady-state solution 
depending on the value of S a . We note in passing that 
a similar phenomenon is known to occur at equilibrium, 
for example at T = in the spin-boson model |241 12"5] . 

The appearance of a bistability is consistent with pre- 
dictions based on a mean field treatment, which is accu- 
rate in the adiabatic limit where ui c —> [I]. For all four 
frequencies studied, the adiabatic effective potentials [26 
shows two distinct minima (upper panel of Fig.|2|, corre- 
sponding to two possible stable configurations |27_. The 
height of the barrier between the two minima is indepen- 
dent of the phonon frequency, however, as clearly evident 
in the figure, the width of the barrier increases as co c de- 



Figure 2: Plots of the occupation of the quantum dot crn(t) 
for different lu c . The solid black, dashed black, solid red, 
and dashed red curves correspond to {n d = 0,S a — 0}, 
{ n d = 1 , S a = 0}, {n d = 0,S a = y / 2ui a \}, and {n d = l,5 a = 
\/2uj a \}, respectively. Upper panel shows the adiabatic ef- 
fective potentials for the different values of uj c . The model 
parameters used are: e d = 0.5eV, and A = 4000cm - . 



creases. This implies that the tunneling time between 
the two configurations also increases as ia c decreases, and 
thus, the bistability (given by the difference between an 
at steady-state for 6 a — and 6 a ^ 0) would increase 
with decreasing uj Ci as indeed is the case. 

The transient dynamics and the approach to steady 
state depends sensitively on the preparation, in partic- 
ular, whether the initial state is close to equilibrium 
({nd — 0, S a — 0} and {rid = 1,6 = \/2uj a X}) or far from 
equilibrium ({n<j = 0,5 a i pha = i/2w a A} and {n d = 1, 
S a = 0}) . For small oj c we observe a rapid decay to steady 
state when the initial preparation is close to equilibrium, 
while for the nonequilibrium preparation, the population 
of the dot <7ii(t) decays to steady state on a time scale 
frT -1 , where T denotes the maximum of leads spectral 
function. As co c increases, the dynamics becomes more 
complex. In particular, a longer time scale in the relax- 
ation of (Tn(t) developes, consistent with the appearance 
of a slow tunneling channel between the two configura- 
tions. However, the existence of this channel does not 
lead to a unique steady state solution for the populations 
at a finite value of ui c , away from the adiabatic limit. 

In Fig. [3] we show <Tn(t) for different values of the 
reorganization energy (polaron shift) A and the dot en- 
ergy Ed- As before, we compare the time dependence of 
<7n(i) for four different initial conditions. In the upper 
panel we show the corresponding adiabatic effective po- 
tentials [26] for the four values of A. For small values of 
A, the bistability clearly disappears (left panels of Fig [3]). 
This is consistent with the adiabatic effective potentials 
showing a single minimum for A < 3000cm -1 (in fact, a 



Figure 3: Similar to Fig.[2]but for different values of A and e d , 
for u) c = 500cm -1 . Upper panel shows the adiabatic effective 
potentials for the different values of A. 



crude estimate based on a mean-field approach suggests 
that below A = 3150cm -1 the bistability vanishes for the 
current parameters). Comparing the relaxation time for 
A = 2000cm" 1 and A = 3000cm" 1 , we find that the latter 
is slower, particularly for the case of S a — y/2uj a X. 

When A is further increased to 4000cm -1 (correspond- 
ing to a nearly symmetric case where the polaron shift 
equals e d ) the relaxation time stretches even more and 
the system decays to a different steady-state depending 
on the value of d a , again consistent with the appearance 
of two stable configurations in the corresponding adia- 
batic effective potential (upper panel of Fig. |3|. While 
the RDM shows a distinct bistability, it is interesting to 
note that this is not the case for the current through the 
quantum dot (not shown), which has a unique steady 
state for the symmetric case (A ~ ej)- 

In the upper right panel of Fig. [3] we show results for 
the case where A = 2000cm -1 and e d = 0.25eV. The 
effective adiabatic potential for this case clearly shows 
two distinct minima, however, the barrier is lower than 
A = 4000cm -1 and e d = 0.5eV. Comparing the two right 
panels of Fig. [3] we find that as A and e d are decreased si- 
multaneously, the bistability decreases and the timescale 
to relax to steady-state also decreases, consistent with 
the adiabatic tunneling picture discussed above. 

In this work, we have combined a RDM formalism with 
the ML-MCTDH-SQR approach to study the nonequi- 
librium dynamics of a many-body quantum system with 
electron-phonon couplings. For a generic model, which is 
widely used to describe phonon-coupled electron trans- 
port in quantum dots and single-molecule junctions, we 
showed that the system may exhibit pronounced bista- 
bility. The analysis reveals that the bistability increases 
for decreasing phonon frequency and depends sensitively 
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on the electron-phonon coupling. This was rationalized 
in terms of the timescales of phonon tunneling between 
the two adiabatic configurations, the phonon relaxation 
time, and the electron relaxation. Based on the RDM 
formalism, we proved that the bistability is associated 
with different initial phonon preparations and not with a 
different initial dot occupations. The present results are 
expected to be of relevance for the interpretation of re- 
cent experiments on charge transport in molecular junc- 
tion, which showed hysteretic behavior [28— 30J , and may 
facilitate further experimental studies in the field of na- 
noclectronics. 
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